%% term_div_Uold_Q
function matrix_A = term_div_Uold_Q(mesh, fespace1, fespace2, Uold)
Nb1 = fespace1.Nb;
Nb2 = fespace2.Nb;
Nlb2 = fespace2.Nlb;
rows = zeros(mesh.Ne*Nlb2,1);
val1 = rows;
flag = 1;
for n = 1:mesh.Ne
    V = mesh.node(mesh.elem(n,:),:);
    [int_points_loc, int_weights_loc] = integral_loc(V, "2d");
    phix = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace1.basis_type, "dx");
    phiy = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace1.basis_type, "dy");
    psi = basis_loc(V, int_points_loc(:,1), int_points_loc(:,2), fespace2.basis_type, "x");
    u1x = phix * Uold(fespace1.elem(n,:));
    u2y = phiy * Uold(fespace1.elem(n,:)+Nb1);
    for i = 1:Nlb2
        f1 = (u1x + u2y) .* psi(:,i);
        r1 = int_weights_loc * f1;
        rows(flag) = fespace2.elem(n,i);
        val1(flag) = r1;
        flag = flag + 1;
    end
end
matrix_A = sparse(rows, 1, val1, Nb2, 1);
end